Solid Earth–atmosphere interaction forces during the 15 January 2022 Tonga eruption

Rapid venting of volcanic material during the 15 January 2022 Tonga eruption generated impulsive downward reaction forces on the Earth of ~2.0 × 1013 N that radiated seismic waves observed throughout the planet, with ~25 s source bursts persisting for ~4.5 hours. The force time history is determined by analysis of teleseismic P waves and Rayleigh waves with periods approximately <50 s, providing insight into the overall volcanic eruption process. The atmospheric acoustic-gravity Lamb wave expanding from the eruption produced broadband ground motions when transiting land, along with driven and conventional tsunami waves. Atmospheric standing acoustic waves near the source produced oscillatory peak forces as large as 4 × 1012 N, exciting resonant solid Earth Rayleigh wave motions at frequencies of 3.7 and 4.6 mHz.

hypothesis that an isotropic implosive source at 5 km depth could account for some of the signal in the displacement waveforms by deconvolving corresponding TGFs with a scalar moment of 10 15 Nm from the data, again using a positivity constraint. The F(t) obtained by assuming an exclusively implosive source can account for initial long-period positive pulse in the P wave arrival for the 0.01-0.05 Hz passband adequately, and it fits the overall signal as well as the single force deconvolution (Fig. S3). The deconvolved moment rate function for the implosive model has general similarity to the point-force time history. This similarity is not unexpected, given the simple shape of the impulse response functions for the two force systems (Fig. S3b,  Fig. S5B), which are essentially flipped in sign. This similarity of the impulse responses leads to ambiguity in the source force system representation, and the deconvolved time series should be viewed as simplified composite time histories of the interfering combination of implosive and reaction forces that generated the P waves.
2) Discussion of tangential component SH waves and Love waves Departures from axisymmetry in the source force system will result in radiation of transverse component SH waves and Love waves. For example, a non-vertical point-force will excite azimuthally varying SH waves with a two-lobed radiation pattern. We examined this possibility by analyzing rotated SH components in the body and surface wave signals. The SH waveforms for most traces had low signal-to-noise ratios (SNR). Traces with good SNR had maximum SH displacements at the arrival time predicted for the USGS-NEIC magnitude 5.8 origin time of <5x10 -7 m, but the data do not robustly resolve a two-lobed radiation pattern. For many stations, SH component ground motion was observed up to 70 s ahead of the expected arrival time, indicating an earlier source. Some of the energy ahead of the theoretical arrival could correspond to signals from a magnitude 4.7 event detected by the USGS-NEIC that occurred 103 s before the magnitude 5.8 origin time. The low SH signal amplitudes and lack of coherent behavior supports assumption of a vertical force model as the primary characterization of the source process for the body waves.
We examined all transverse component recordings for long period Love wave motions as well. Again, transverse component motions are apparent, with amplitudes that are up to about 75% as large as the Rayleigh wave signals. The Love wave observations can be accounted for by up to 30° deflection of the point-force from vertical, directed toward the WNW, but the dominant period range of the observations is from 10 to 50 s, so deflections from the great circle path and similarity of the group velocities for these periods may cause contamination of the putative transverse components. Modest deflection of a point-force from the vertical produces minor azimuthal pattern in the P and Rayleigh wave radiation, and the amplitude data for those phases do not confirm a geometry that accounts for the tangential motions, so we focus on axisymmetric vertical point force and implosive sources for the source representation.

3) Model dependence of F(t): AK135, PREM_a
We explored the robustness of our surface waves results by removing the point force impulse response for two different standard Earth models, AK135 and PREM_a. The latter model considers an anisotropic Earth structure and has slightly sharper changes of velocity than AK135 in the first 600 km. Fig. S8 shows examples of vertical component (Rayleigh wave) impulse responses for AK135 and PREM_a and corresponding signals filtered in the 0.01-0.05 Hz passband (used for the deconvolutions) and the 0.01-0.03 Hz passband. In the filtered passbands the strong differences in short-period Rayleigh wave dispersion for the two models are greatly reduced, however, there are still minor differences in the Green's function waveforms, which account for the differences in the deconvolutions. The main effect is that inaccuracies in the predicted dispersion across the dominant passband of the data from 0.025-0.05 Hz result in sidelobes to primary pulses associated with the well-predicted source radiation. The sidelobes can be skewed early or late depending on the trend in dispersion discrepancy.
The median stacks for the AK135 and PREM_a models, not shifted to absolute time, are shown in Fig. S9 and Fig. S10, respectively. The comparison of ~16,500 s long F(t) shows that we can confidently resolve broad features of the eruption using the two models (Figs. S9A, S10A). The first stage of the eruption is composed of a series of force pulses that decay in amplitude with time lasting ~5000 s. The most vigorous reaction forces occurred during the first ~700 s. A secondary, and much less energetic, stage occurred at ~15,000 s and it sustained for ~1,000 s, making the total time of the eruption process ~16,000 s long.
Finer scale details of the first 700 s of the eruption ( Fig. S9B and Fig. S10B) show that there are some differences in the stacked F(t) for the two Earth models. The Rayleigh wave energy is relatively narrow band, with most power in the period band 20 to 40 s. This relatively shortperiod band is strongly influenced by the crustal and uppermost mantle structure and the dispersion curves at these periods differ for models AK135 and PREM_a. No one-dimensional Earth model is expected to predict the dispersion accurately for the various paths sampled by the data, and errors in predicted dispersion in the impulse responses for the narrowband data will lead to side-lobes in the deconvolutions. This appears to cause some early artifacts in the PREM_a deconvolutions and likely some later artifacts in the AK135 deconvolutions, relative to common main pulses. Nevertheless, the deconvolutions for the two models show the same main characteristics of the force time history as obtained from the P waves deconvolution for the early part of the eruption. The peak force is followed by a long-period pulse, which is divided into two different pulses in the AK135 surface wave stack. The second set of three force pulses ~100s after the long-period pulse observed in the P wave deconvolution can also be seen in the two surface waves stacks, but the PREM_a F(t) seems to have some early side lobes. The detailed velocity structure of the real Earth is not known sufficiently to account for the short-period surface wave dispersion precisely, so the estimates of F(t) complexity are not uniquely resolved. The very different structures of AK135 and PREM_a likely span the range of structural mismatches expected for 1D models, so we have confidence in the common features of the deconvolutions. It appears that the main characteristics of the force system acting on the volcano as the eruption develops are captured quite well by the long surface wavetrains radiated from the source.
Fig. S1. -Corrected P amplitudes for the first large positive ground displacement pulse versus azimuth. A. Peak amplitudes of the instrument-corrected ground displacement signals of the January 15, 2022 Tonga eruption in the passband 0.01-0.05 Hz for the first large P wave arrival at global broadband seismic stations in the epicentral distance range 50° to 90° (black dots), with corrections for geometric spreading to a distance of 78.5° (orange dots), show little azimuthal pattern. The median values of the observations for different passbands are used to scale the normalized median stacks of bin-averaged P wave ground displacements (Fig. S2) giving the true-amplitude median ground motion stacks shown in Figs. S3 and S4. B. Peak amplitudes of the instrument-corrected ground displacement signals of the November 11, 2022 MW 7.3 earthquake in the Tonga Islands in the passband 0.01-0.05 Hz for the first large P wave arrival at global broadband seismic stations in the epicentral distance range 50° to 90°, with corrections for geometric spreading and radiation pattern to a distance of 78.3°. The radiation pattern and geometric spread corrected amplitudes for the MW 7.3 earthquake also show little azimuthal variation, with similarities with the eruption data indicating that the residual pattern is primarily due to attenuation and receiver effects and not to source effects.

Fig. S5 -Deconvolution of the P wave stack by an implosion impulse response. A:
Global stack for P wave ground displacement filtered in the 0.01-0.05 Hz band at 78.5° epicentral distance (same as Fig. S3A). The red dotted line indicates the expected P wave arrival time for the magnitude 5.8 event reported by the U.S. Geological Survey National Earthquake Information Center, as in Fig. 3. B: Synthetic P wave impulse response seismogram at 78.5° epicentral distance for an isotropic implosive point-source with seismic moment of 10 15 N-m located at 5 km depth in model AK135, filtered in the 0.01-0.05 Hz passband. C: Moment rate source-time function obtained by iterative deconvolution with positivity constraint. D: Comparison between the observed P waveform (black curve) and the reconstituted waveform (red curve) obtained by convolving the implosion moment rate function time series (C) with the impulse response seismogram (B).    Individual station deconvolutions are shown with color traces with the legend on the right. The black curve is the median stack of the aligned F(t) estimates. B: The initial 700s of the aligned F(t), with individual stations shown by color lines and the median stack shown by the black line. Note that different stations were screened for stacking F(t) for different time intervals due to differences in signal-to-noise ratios during the ~4.5 hr long eruption.The traces are not shifted to absolute time and the timing here depends on the relative alignment of the F(t) obtained with the AK135 model. Individual station deconvolutions are shown with color traces with the legend on the right. The black curve is the median stack of the aligned F(t) estimates. B: The initial 700s of the aligned F(t), with individual stations shown with color lines and the median stack shown by the black line. Note that different stations were screened for stacking F(t) for different time intervals due to differences in signal-to-noise ratios during the ~4.5 hr long eruption. The PREM_a Earth model causes some side-lobes in the deconvolution due to slight errors in Rayleigh wave dispersion in the data passband, but the major features are similar to Fig. S9. The traces are not shifted to absolute time and the timing here depends on the relative alignment of the F(t) obtained with the PREM_a model.